



GIS714: geospatial computation and simulation
Graduate level course required for Geospatial Analytics PhD students
Previously taught with GRASS tutorials, instructions posted on webpage
Goal: convert to GRASS-based assignments in Jupyter Notebooks


“The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text.” – jupyter.org

The default kernel is IPython.
import sys
v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")
We are using Python 3.9.5
Output can be interactive.
from ipywidgets import interact
def f(x):
return x
interact(f, x=10);
...enhances the existing GRASS Python API to allow Jupyter Notebook users to easily manage GRASS data, visualize data including spatio-temporal datasets and 3D visualizations
GIS714: Geospatial Computation and Simulation, Spring 2022
We'll use an example from the course on water simulation...
# Import Python standard library and IPython packages we need.
import subprocess
import sys
# Ask GRASS GIS where its Python packages are.
sys.path.append(
subprocess.check_output(["grass80", "--config", "python_path"], text=True, shell=True).strip()
)
# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj
# Start GRASS Session
session = gj.init("data", "nc_spm_08_grass7", "PERMANENT")
# Make a new mapset for this assignment
gs.run_command("g.mapset", mapset="HW3_water_simulation", location="nc_spm_08_grass7", flags="c")
# Set region to the high resolution study area
gs.run_command("g.region", region="rural_1m")
And then we are ready to begin using GRASS!
Since grass.jupyter is written in Python, we'll use the GRASS Python API.
We create a map of flooded area with r.lake (GRASS manual for r.lake) by providing a water level and a seed point:
gs.run_command("r.lake", elevation="elev_lid792_1m", water_level=113.5, lake="flood1", coordinates="638728,220278")
# See results
flood1 = gj.Map(use_region=True, width=200)
flood1.d_rast(map="elev_lid792_1m")
flood1.d_rast(map="flood1")
# Display map
flood1.show()
Increase water level to 113.7m and 114.0m and create flooded area maps at these two levels.
#### Your Answer Here
We can also make 3D maps using similar syntax
elev_map = gj.Map3D()
elev_map.render(elevation_map="elevation", color_map="elevation", perspective=20)
elev_map.overlay.d_legend(raster="elevation", at=(60, 97, 87, 92))
elev_map.show()
Calling display modules using Python magic short-cut:
d.rast -> Map.d_rast(), Map3D.d_rast
We will use two GRASS addons, r.stream.distance and r.lake.series, to estimate inundation with Height Above Nearest Drainage methodology (A.D. Nobre, 2011).
For this section, we will change our computation region to elevation which is a larger study area than we used above. Because we are using a new region and a higher threshold of 100,000, we need to run r.watershed again to compute the flow accumulation, drainage and streams.
gs.run_command("g.region", raster="elevation")
gs.run_command("r.watershed", elevation="elevation", accumulation="flowacc", drainage="drainage", stream="streams_100k", threshold=100000)
gs.run_command("r.to.vect", input="streams_100k", output="streams_100k", type="line")
Now we use r.stream.distance with output parameter difference to compute new raster where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains. This is our HAND terrain model.
gs.run_command("r.stream.distance", stream_rast="streams_100k", direction="drainage", elevation="elevation", method="downstream", difference="above_stream")
With r.lake.series, we can create a series of inundation maps with rising water levels. r.lake.series creates a space-time dataset. We can use temporal modules to further work with the data...
gs.run_command("r.lake.series", elevation="above_stream", start_water_level=0, end_water_level=5, water_level_step=0.5,
output="inundation", seed_raster="streams_100k")
... or visualize the flood with TimeSeriesMap.
flood_map = gj.TimeSeriesMap(use_region=True)
flood_map.d_rast(map="elevation")
flood_map.add_raster_series("inundation")
flood_map.d_legend(color="black", at=(2,30,2,5)) #Add legend
flood_map.show()
HBox(children=(Play(value=0, description='Press play', interval=500, max=10), SelectionSlider(description='Dat…
Interactivity is often useful for explore geospatial data.
folium...

Moving data from GRASS GIS location projection to folium is a challenge
folium is projected in Web Mercator (EPSG:3857)
However, any coordinates (i.e. vector data) need to be specified in degrees of latitude and longitude (WGS84, EPSG:4326)
folium reprojects latitude and longitude to Web Mercator internally
We pass data to folium by reprojecting to temporary locations

import folium
# Create figure
fig = folium.Figure(width=600, height=400)
# Create a map to add to the figure later
m = folium.Map(tiles="Stamen Terrain", location=[35.761168,-78.668271], zoom_start=13)
# Create and add elevation layer to map
gj.Raster("elevation", opacity=0.5).add_to(m)
# Do some cool folium stuff!
# Like make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
[35.781608,-78.675800], popup="<i>Center For Geospatial Analytics</i>", tooltip=tooltip
).add_to(m)
# and a circle
folium.Circle(
radius=120,
location=[35.769781,-78.663160],
popup="Great Picnic Area",
color="crimson",
fill=False,
tooltip=tooltip
).add_to(m)
# Add the map to the figure
fig.add_child(m)
# Display figure
fig
The InteractiveMap class allows users unfamiliar with folium to produce maps easily.
# Create Interactive Map
fig = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
fig.add_raster("elevation", title="Elevation Raster", opacity=0.8)
#fig.add_vector("roadsmajor")
fig.add_layer_control(position = "bottomright")
# Display map
fig.show()
Notebook format was generally well received:
We did experience some challenges:
FOSS4G '22 GRASS GIS Workshop Materials (Anna Petrasova)
GIS714: Advanced Geospatial Computation and Simulation Course Repository
Helena Mitasova
Stefan Blumentrath
Vero Andreo
GIS714 Class Spring 2021
… and the GRASS community, NC State Center for Geospatial Analytics, Google Summer of Code